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INTRODUCTION 


Boundary  element  methods  (BEM)  have  become  an  accepted  alternative 
to  domain-based  methods  (finite  elements  and  finite  difference)  for  many 
classes  of  boundary  value  problems  (Refs  1  and  2).  The  direct  (DBEM) 
and  indirect  (IBEM)  formulations  are  two  well  established  formulations 
that  are  used  to  solve  a  variety  of  boundary  value  problems  in  engineer¬ 
ing.  The  key  to  these  formulations  is  the  existence  of  a  fundamental 
singular  solution  (i.e.,  the  solution  for  a  point  source  or  point  load 
in  an  infinite  domain  --  a  free-space  Green's  function). 

This  singular  solution  is  jointly  the  blessing  and  curse  of  the 
methods.  Since  it  is  a  solution  of  the  governing  differential  equations, 
problems  defined  by  a  self-adjoint  differential  operator  can  be  expressed 
by  integral  equations.  Numerically  this  means  that  for  many  classes  of 
problems  BEMs  require  only  the  boundary  to  be  discretized.  The  dimension 
ality  of  problems  is  reduced  by  one,  and  thus  the  methods  lend  themselves 
to  modeling  infinite  domains.  The  reduction  in  dimensionality  reduces 
the  number  of  degrees  of  freedom  in  the  approximation;  however,  the  re¬ 
sulting  equations  are  nonsymmetric  and  fully  populated  in  contrast  to 
the  symmetric  sparse  systems  characteristic  of  domain-based  methods. 

The  singular  nature  of  the  fundamental  solution  allows  the  methods  to 
effectively  approximate  problems  having  singular  or  high-gradient  stress 
fields.  This  same  singular  nature  complicates  the  numerical  integrations 
inherent  to  the  methods.  The  fundamental  singular  solution  is  well  known 
for  many  classes  of  boundary  value  problems  with  homogeneous  domains. 
However,  for  an  arbitrary- inhomogeneous  domain  the  fundamental  solution 
is  not  known,  and  additional  numerical  effort  is  required  (Ref  3). 

Since  the  methods  are  relatively  new,  they  have  not  been  as  extensively 
developed  as  domain  methods  for  nonlinear  applications.  They  also  lack 
the  generality  (in  terms  of  extensive  continuum  and  structural  element 
librai ies)  that  commercial  finite  element  computer  programs  possess. 


The  finite  element  method  (FEM)  has  been  the  most  highly  developed 
approximate  solution  technique  for  many  classes  boundary  value  problems 
(Ref  39).  The  most  general  derivation  follows  from  a  weighted  residual 
approximation  of  the  boundary  value  problem.  The  error  introduced  by  an 
approximate  solution  is  orthogonalized  with  respect  to  a  set  of  weighting 
functions  to  obtain  the  "best"  approximate  solution.  For  the  Galerkin 
weighted  residual  statement  the  weighting  functions  correspond  to  the 
basis  functions  in  which  the  approximate  solution  is  expanded.  For  prob¬ 
lems  with  a  self-adjoint  operator,  the  weak  form  of  the  Galerkin  statement 
is  equivalent  to  the  variational  statement  of  the  problem. 

Basis  functions  in  the  FEM  usually  are  locally  based  polynomials 
for  which  the  integrations  are  relatively  simple.  The  local  nature  of 
these  functions  leads  to  a  sparse  system  of  equations,  but  this  same 
attribute  hinders  its  ability  to  model  infinite  domains.  To  compensate 
for  this  limitation,  special  elements  (infinite  elements)  have  been  de¬ 
veloped  to  model  infinite  domains  which  include  a  shape  function  based 
on  the  form  of  analytical  solutions  (Ref  4).  This  contrasts  with  the 
BEM,  where  rather  than  using  analytical  solutions  to  determine  the  form 
of  shape  functions,  the  analytical  solutions  are  used  directly  as  weighting 
functions  (Ref  5).  The  generality  of  the  Galerkin  and  variational 
statements  allows  nonlinearities  and  domain  inhomogenet ies  to  be 
accommodated.  For  most  structural  applications  the  resulting  system  of 
equations  is  also  symmetrical. 

Considering  the  individual  strengths  of  boundary-  and  domain-based 
methods,  some  classes  of  problems  may  be  most  effectively  solved  by  com¬ 
bining  the  methods.  Nonlinear  soil-structure  interaction  and  fracture 
mechanics  are  two  applications  where  a  combined  solution  approach  can 
potentially  provide  a  more  effective  analysis.  For  nonlinear  soil- 
structure  interaction,  the  FEM  could  be  used  to  model  the  structure  and 
a  region  of  the  soil  which  is  expected  to  exhibit  nonlinear  behavior; 
the  BEM  could  be  used  to  approximate  the  infinite  domain.  For  fracture 
mechanics,  the  FEM  might  be  used  to  model  all  of  the  problem  except  a 
region  local  to  the  crack  tip  where  the  BEM  could  provide  a  crack  tip 
"element"  capable  of  representing  very  high  gradients. 
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Objective 


The  ultimate  objective  of  this  research  is  to  determine  whether  a 
combination  of  boundary  element  and  finite  element  solution  methods 
could  more  effectively  solve  nonlinear  structural/geotechnical  problems 
than  single  method  approaches.  The  immediate  emphasis  is  on  the  formula¬ 
tion  of  stiffness  matrices  for  both  the  indirect  and  direct  boundary 
element  methods  (IBEM  and  DBEM,  respectively). 

Background 

The  idea  of  coupling  the  FEM  and  BEM  is  attributed  to  Wexler  (Ref  6), 
who  used  integral  equation  solutions  to  constrain  FE  solutions  for  infinite 
domain  field  problems  in  the  early  1970s.  Other  early  examples  of  using 
a  combined  solution  approach  dealt  with  wave  phenomena,  Chen  and  Mei 
(Ref  7)  and  Shaw  (Ref  8).  Zienkiewicz  et  al.  (Ref  9)  proposed  a  combined 
solution  approach  for  statics  problems  followed  by  Osias'  use  of  a  coupled 
solution  to  solve  elastostatic  problems  (Ref  10). 

There  are  numerous  approaches  to  coupling  the  methods  (Ref  11).  In 
this  study  the  main  classification  is  based  upon  "which  method  hosts  the 
coupling."  For  a  BEM-hosted  coupling  the  FEM  subdomain  is  treated  as  a 
BE  region;  equilibrium  and  compatibility  are  approximately  enforced  con¬ 
tinuously  along  the  interface,  analogous  to  the  BEM  approach  to  modeling 
piece-wise  homogeneous  problems.  The  resulting  equations  are  nonsymmet- 
ric,  and  thus  this  approach  only  lends  itself  to  problems  numerically 
dominated  by  the  BE  subdomain  (Refs  11,  12,  13,  14  and  15).  For  a  FEM- 
hosted  coupling  the  BEM  subdomain  is  treated  as  one  or  many  finite  ele¬ 
ments.  In  this  "super-element"  approach  the  resulting  BE  equations  re¬ 
flect  the  "equivalent"  stiffness  of  the  subdomain  and  can  be  directly 
assembled  by  a  FE  program.  The  term  "super-element"  reflects  the  high 
degree  of  connectivity  which  will  typically  be  present  for  a  BE  region. 

A  coupled  solution  approach  does  not  seem  to  necessitate  a  displacement- 
based  FE  formulation;  however  due  to  its  commercial  success,  it  appears 
all  the  studies  of  coupled  solutions  have  been  based  on  this  formulation. 
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Within  FEM-hosted  couplings  the  next  level  of  classification  is  the 
procedure  used  to  obtain  the  stiffness  relationship.  The  two  approaches 
are:  (1)  a  variational  approach  where  a  boundary  variational  equation  is 
the  basis  of  the  relationship,  and  (2)  a  direct  approach  where  the  BE 
equations  are  manipulated  into  a  stiffness  form.  A  key  step  to  both 
approaches,  using  either  BE  formulation,  is  establishing  a  relationship 
between  the  nodal  tractions  and  the  nodal  displacements;  an  exception  to 
this  is  a  variational  approach  with  the  IBEM  which  retains  the  artificial 
tractions  as  unknowns  in  the  system  of  equations  (Refs  11  and  12). 

Mei  (Ref  16)  was  apparently  the  first  to  formulate  a  coupled  solution 
based  on  variational  techniques  (Ref  13).  Other  early  efforts  using  a 
variational  approach  include:  Jeng  and  Wexler  (Ref  17);  Shaw  (Ref  18); 
Zienkiewicz,  Kelly,  and  Bettess  (Ref  9);  and  Kelly,  Mustoe,  and  Zienkiewicz 
(Ref  11).  Stationarity  of  the  variational  statement  necessarily  produces 
a  symmetric  system,  and  thus  this  approach  is  sometimes  grouped  with 
some  direct  approaches  as  "symmetric  couplings."  In  the  variational 
approach,  a  variational  statement  of  the  problem  is  written  in  a  boundary 
form  with  the  absence  of  body  sources  or  forces  (Refs  11  and  19).  For 
elastostatics  the  variational  statement  corresponds  to  a  boundary  form 
of  the  potential  energy  functional.  There  have  been  many  modifications 
of  this  functional  for  various  purposes  including  (1)  to  relax  essential 
boundary  conditions  or  interface  compatibility  and  (2)  to  enforce  equilib¬ 
rium  (Ref  19  and  20). 

With  a  functional,  relationships  between  boundary  values  obtained 
from  the  integral  equations  are  substituted  into  the  functional,  reduc¬ 
ing  the  number  of  unknowns.  For  a  displacement-based  FE,  formulation 
the  nodal  tractions  are  usually  expressed  in  terms  of  the  nodal  displace¬ 
ments,  which  requires  the  inversion  of  a  fully  populated  coefficient 
matrix.  The  previously  mentioned  exception  to  this  approach  (using  the 
TBEM),  writes  the  boundary  tractions  in  terms  of  the  artificial  tractions. 
In  this  case  stationarity  is  taken  with  respect  to  both  nodal  displace¬ 
ments  and  artificial  boundary  tractions  instead  of  nodal  displacements 
alone.  The  lack  of  a  matrix  inversion  is  paid  for  by  an  increase  in  the 
number  of  unknowns  and  a  double  integration  process. 
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In  contrast  to  the  variational  formulations,  Brebbia  (see  e.g., 

Ref  2,  5,  or  21)  directly  manipulated  the  BE  equations  into  a  stiffness 
form.  The  Maxwell-Betti  reciprocal  theorem  proves  that  a  stiffness 
matrix  should  be  symmetric  for  a  linear  elastic  domain;  however,  the 
direct  manipulation  does  not  produce  a  symmetric  system.  Brebbia 
symmetrized  the  system  by  averaging  the  off-diagonal  terms  based  on  an 
error  minimization  argument  using  the  method  of  least  squares.  The 
resulting  stiffness  matrix  was  of  the  same  form  as  obtained  by 
variational  methods  using  the  potential  energy  functional  (Ref  11). 

Despite  the  apparent  agreement  between  the  two  approaches,  agreement 
between  proponents  of  the  approaches  is  scarce.  Kelly  et  al.  (Ref  11) 
state  that  the  modification  to  obtain  symmetry  "apparently  has  no  valid 
foundation  except  that  the  energy  approach  led  to  a  similar  form."  Brebbia 
et  al.  (Ref  2)  later  admitted  that  there  is  no  rationale  to  justify  the 
error  minimization  and  that  the  matrix  should  in  fact  be  nonsymmetr ic . 

In  addition,  he  questions  the  combination  of  variational  principles  and 


integral  equations  used  in  the  variational  approach.  His  argument  is  not 
implying  that  a  nonsymmetr ic  stiffness  matrix  is  physically  meaningful 
for  an  elastic  domain;  rather  he  seems  to  be  indicating  that  the  numerical 
approximations  used  by  the  BEM  will  inherently  produce  a  nonsymmetric 
stiffness  matrix.  Initial  efforts  indicated  that  the  symmetric  matrix 
provided  better  results;  however,  more  recent  comparisons  (Ref  20)  using 
quadratic  elements  and  providing  for  geometric  discontinuities  indicate 
the  contrary.  Fortunately  for  cases  which  have  a  relatively  fine  boundary 
discretization,  Roudolphi  (Ref  22)  has  found  the  degree  of  nonsymmetry 
to  be  insignificant.  This  is  likely  to  be  the  case  for  geotechnical 
app 1 icat ions . 

Another  problem  which  is  common  to  both  approaches  is  failure  of 
the  stiffness  matrix  to  satisfy  equilibrium.  Equilibrium  requires  the 
terms  in  each  column  corresponding  to  a  given  global  direction  to  sum  to 
zero.  With  the  DBEM  the  system  of  equations  can  be  augmented  by  two 
additional  equations  which  cause  the  approximate  traction  distribution 
to  identically  satisfy  equilibrium  (Ref  11).  By  this  approach  the 
equilibrium  problem  is  addressed  while  establishing  the  nodal  traction- 
displacement  relationship  instead  of  after  the  stiffness  matrix  is 
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formed.  Hartman  (Ref  23)  gives  a  theoretical  discussion  on  the  symmetry 
and  equilibrium  problems  and  provides  advice  on  how  to  correct  them. 
Tullberg  and  Belteus  (Ref  20)  compare  the  accuracy  and  convergence 
characteristics  of  seven  stiffness  formulations  for  two  simple 
two-dimensional  problems  --  uniaxial  tension  and  bending.  The  seven 
stiffness  matrices  were  obtained  using  the  DBEM  and  included  all  the 
variations  mentioned  above. 

In  addition  to  the  differences  used  here  to  classify  formulations, 
there  have  been  numerous  other  variations.  Some  researchers  (Ref  19) 
have  used  element  condensation  to  reduce  the  unknowns  for  the  BE  super¬ 
element  to  those  located  on  the  interface  between  the  two  methods.  This 
could  be  useful  in  fracture  mechanics  applications  where  numerous  boundary 
elements  might  be  used  near  the  crack  tip.  As  previously  mentioned, 
modified  functionals  have  been  use  1  to  relax  compatibility  requirements 
(e.g.,  when  interfacing  subdomains  with  different  shape  functions). 

Variations  in  the  BE  formulations  themselves  introduce  additional 
classification  parameters.  Though  not  addressed  in  this  study,  for  some 
applications  traction  discontinuity  is  an  important  consideration.  Many 
ways  of  treating  this  problem  have  been  reported  (Ref  20,  25,  26,  and 
27).  A  simple  though  less  rigorous  approach  is  to  position  the  colloca¬ 
tion  points  at  the  interior  of  the  element.  When  all  of  the  collocation 
points  associated  with  an  element  are  internally  positioned,  the  elements 
are  referred  to  as  "discontinuous"  or  "nonconforming"  elements  (Ref  28). 
Though  the  discontinuous  elements  add  some  flexibility  to  the  method  in 
terms  of  mesh  refinement,  it  is  at  the  expense  of  increasing  the  number 
of  equations  by  up  to  100  percent.  Discontinuous  elements  would  appear 
to  complicate  a  coupling  algorithm  since  the  collocation  points  do  not 
correspond  to  node  points. 

Relationship  to  Previous  NGEL  Work 

The  primary  role  of  this  research  of  the  BEM  has  been  to  compliment 
FEM  analysis  capabilities  (i.e.,  coupled  solution  techniques).  The  BEM 
is  in  general  a  more  analytically  rigorous  technique,  and  thus  the  fund¬ 
amentals  of  the  method  have  been  investigated.  As  a  result  of  these 
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efforts  the  potential  of  the  method  as  a  separate  analysis  tool  has  also 
been  observed.  Thus  the  efforts  at  coupling  have  emphasized  a  modular 
approach  where  it  is  assumed  that  both  analysis  methods  can  serve  as  an 
individual  or  a  coupled  tool. 

In  NCEL's  first  study  (Ref  29)  the  IBEM  formulation  was 
investigated  in  one-  and  two-dimensional  elastostatics  and  numerical 
results  were  compared  with  the  DBEM  and  analytical  solutions.  A  crude 
element  (constant  distribution)  was  used  but  the  potential  of  the  method 
was  evident.  This  study  indicated  the  need  for  higher  order  elements 
and  integration  techniques  to  improve  results  in  the  near  boundary 
region. 

In  NCEL's  next  study  (Ref  15)  a  coupling  of  the  BE  and  FE  methods 
was  investigated.  A  BEM-hosted  coupling  was  formulated,  approximating 
compatibility  and  equilibrium  on  the  interfaces.  Qualitatively  the 
coupling  was  a  success,  but  the  simple  constant  elements  used  for  both 
methods  prevented  substantial  quantitative  evaluation.  For  nonlinear 
soil-structure  applications  where  the  FEM  is  used  to  model  the  nonlinear 
region,  the  FEM  should  host  the  coupling.  Two  areas  required  further 
study:  (1)  accuracy  of  the  BEM  formulations,  and  (2)  a  FEM-hosted 
coupling  approach. 

In  the  third  study  (Refs  30  and  31)  higher  order  isoparametric  bound¬ 
ary  elements  were  investigated  along  with  their  associated  integration 
techniques.  Lachat  and  Watson  (Ref  32)  adapted  the  isoparametric  element 
formulation  directly  from  FE  technology.  However,  for  the  BEM  the  inte¬ 
grands  in  the  coefficient  calculations  are  singular  and  thus  their  integra¬ 
tion  requires  special  attention.  Though  much  effort  had  been  devoted  to 
dealing  with  the  integration  over  the  singularity  an  effective  numerical 
approach  to  performing  the  integrations  in  the  near  boundary  region  was 
lacking.  A  recursive  algorithm  which  adaptively  subdivides  the  element 
to  maintain  consistent  accuracy  independent  of  the  response  point  location 
was  developed.  The  new  element  formulation  has  effectively  eliminated 
near  boundary  error  resulting  from  coarse  numerical  integration  and  has 
provided  valuable  insight  to  other  errors  inherent  to  the  boundary 
element  methods. 
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Scope 


In  this  study,  different  approaches  to  obtaining  a  FEM-hosted  coupling 
are  investigated.  The.  emphasis  is  on:  (1)  conceptually  unifying  the  dif¬ 
ferent  coupling  techniques  and  (2)  providing  an  algorithm  to  facilitate 
future  implementation.  This  preliminary  report  on  the  study  emphasizes 
the  development  of  an  algorithm  for  obtaining  a  stiffness  matrix  from 
RF.Ms.  A  future  report  should  provide  a  more  detailed  theoretical  back¬ 
ground  and  numerical  results.  A  physically  intuitive  development  of  an 
I  REM  stiffness  matrix  is  included,  but  the  DBEM  is  also  addressed.  Over¬ 
views  of  both  the  DBEM  and  IBEM  formulations  are  given  in  the  following 
sections  to  establish  nomenclature  and  to  present  the  necessary  equations 
for  background. 

Nomenc 1 ature 

A  list  of  principal  symbols  used  in  the  report  is  included  as 
Appendix  A.  Additional  variables  which  have  a  limited  scope  of  use  are 
defined  within  the  text.  Modifications  (such  as  the  addition  of  a  sub¬ 
script)  are  also  defined  in  the  text  unless  listed  in  Appendix  A. 

Instead  of  the  matrix  notation  defined,  indicial  notation  will  sometimes 
be  used  to  clarify  equations.  Indices  will  always  be  lower  case  while 
identifying  subscripts  or  superscripts  will  be  upper  case. 


THEORY 


1'his  study  focuses  on  the  development  of  a  stiffness  matrix  for  a 
coupled  solution  approach.  Each  method  of  analysis  is  applied  to 
port,  ions  of  the  domain  where  it.  is  best  suited.  For  illustration  the 
domain  Q  ( Figure  1)  is  divided  into  two  subdomains  51^  ,  the  FF.  subdomain, 
and  ,  the  RE  subdomain.  The  equations  of  equilibrium  are  given  by 
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and  the  compatibility  equations  are  given  by 


t  .  .  =  r(u.  .  +  u .  .) 

ij  2  i,j  J,i 


The  BE  subdomain  is  assumed  to  be  a  linear  isotopic  homogeneous  elastic 
material.  Thus  the  governing  constitutive  relations,  generalized  Hook's 
law,  are  given  as 


a  .  .  =  2  y  £  .  .  +  X  e,  .  6  .  . 
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where  X  and  y  are  Lame's  constants  expressed  in  terms  of  Young's  modulus 
(E)  and  Poisson's  ratio  (v)  as 


(1  +  v)(l  -  2  v)  ’ 


2(1  +  v) 


Alternatively  the  equilibrium  (1),  compatibility  (2),  and  constitutive 
(3)  relations  can  be  combined  to  give  a  single  set  of  equations  in  terms 
of  displacement  written  as 


U  u  .  ,  .  +  f  X  4-  y )  u  .  .  .  4-  ip  .  =  0 
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the  Navier  equations.  The  boundary  conditions  are  given  by 


u.(x)  =  u  (x)  on 


t,(x)  =  t.(x)  on  F,r 


where  u,(x)  and  t^(x)  are  prescribed  distributions  of  boundary 
displacements  and  tractions,  respectively.  The  simple  notation  does  rot 
imply,  that  the  boundary  conditions  are  mutually  exclusive;  the  fully 
mixed  boundary  value  problem  is  addressed. 
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The  key  solution  to  Equation  (5)  for  the  more  popular  BEMs  is  the 
fundanental  singular  solution.  This  is  the  solution  due  to  a  concen¬ 
trated  load  in  an  infinite  space  which  has  the  same  dimension  as  the 
problem  to  be  solved.  For  plane  strain  conditions  the  Kelvin  solution 
expresses  the  displacement  field  u^(x)  due  to  a  unit  force  e^CU  in  an 
infinite  plane.  The  indices  i,  j,  and  k  assume  values  of  1  or  2,  and 
repeated  indices  imply  summation.  The  Kelvin  solution  is  given  by 


ui(?)  =  Gik(*>4)  ek(P 
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where 
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arbitrary  constant  tensor  based  on  zero  displacement 
reference  distance 
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By  incorporating  Equations  (2)  and  (3),  the  stress  field  o„(x)  is  given 


as 


°ii(-}  =  Tiik(?’^ 
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where 


Tuk'?'5> 


c,  (6  y  .  +  6  ..  y  .  -  6  .  . y 
4  ik' j  jk  r  ly  k 


i  +  2yiyjyk 
.•  ;yJ  + - o 


(8b) 


'3  4n( 1  -  v) 


2v 


10 


3 


Equilibrium  conditions  applied  at  a  boundary  point,  with  a  unit 
outward  normal  n^Cx),  and  Equation  (8a)  combine  to  give  the  surface 


tractions  t^(&)  as 


V*)  =  FlkCx,4)  e^(^)  (9a) 

where 

Flk(?.5>  =  4  VVi  -  vk)  +  (  C  46  lk+  )  V)  (9b) 

Figures  2  and  3  illustrate  the  singular  behavior  of  the  fundamental 
solution  and  T^  respectively,  for  a  point  load  applied  at  the 
origin  of  the  Cartesian  system.  The  singularity  of  G  is  order  ln(r), 
while  the  singularity  of  T,  obtained  from  derivatives  of  G,  is  order 
1/r.  The  plane  strain  solution  can  be  converted  to  the  plane  stress 
solution  by  specifying  an  effective  Poisson  ratio  v  =  v/(l+v). 

The  fundamental  solution  is  a  key  ingredient  in  formulating 
integral  equations  for  the  indirect  and  direct  BEMs.  Integral  equations 
are  an  equivalent  statement  of  a  boundary  value  problem.  Finite 
difference  solutions  approximate  the  differential  equations;  finite 
element  solutions  approximate  the  stationarity  of  the  variational 
statement  (in  some  instances);  and  BEMs  approximate  the  integral 
equations . 

The  following  subsections  provide  an  overview  of  both  the  direct 
and  indirect  BEMs  and  their  associated  integral  equations.  These  over¬ 
views  are  followed  by  a  derivation  of  an  IBEM  stiffness  matrix.  The 
last  subsection  discusses  the  discontinuous  nature  of  coupled  solution 
techniques--an  inherent  trait  of  integral  equation  methods. 

Overview  of  Direct  Boundary  Element  Method 

The  direct  boundary  element  method  is  the  most  highly  developed  of 
all  the  integral  equation  methods.  It  was  first  applied  to  elastostatics 
by  Rizzo  (Ref  33).  Cruse  and  Rizzo  (Ref  3U)  followed  with  a  solution  of 


V  A  .v .VA  AAA  ,■>  ,V; 
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the  general  transient  el asto-dynamic  problem.  Since  these  early  applica¬ 
tions,  it  has  matured  numerically  and  the  scope  of  application  has 
significantly  broadened. 

The  integral  equations  on  which  the  DBEM  is  based  are  known  as  the 
Somigliana  identities.  The  first  identity  has  been  obtained  by  two  dif¬ 
ferent  derivations.  In  both  cases  the  fundamental  solution  plays  a  key 
role.  The  first  derivation  uses  Betti's  reciprocal  work  theorem  in  which 
one  system  is  the  actual  boundary  value  problem  and  the  second  system 
corresponds  to  the  fundamental  solution  (Ref  33).  Betti  (1872-73)  and 
Somigliana  (1885-86)  were  the  first  to  apply  potential  methods  to  elas¬ 
ticity  (Ref  19).  Thus  it  is  not  surprising  that  Betti's  theorem  and 
Somigliana' s  first  identity  are  the  Navier  equation  counterparts  to 
Green's  second  and  third  formula  from  potential  theory. 

Somigliana' s  first  identity  is  given  by 


U-(P  =  ^(x)  ‘  Fi^(x,|)  u.(x)]  dx 
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+  | G  (z,§)  *A(z)  dz  (10a) 
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in  which  the  displacement  field  is  written  in  terms  of  the  boundary  tractions, 
boundary  displacements,  and  body  forces.  The  second  identity  is  obtained 
by  combining  Equation  (10a)  with  the  compatibility  (Equation  (2))  and 
the  constitutive  (Equation  (3))  relations  giving  the  stress  tensor  as 

ajk(^  =  jfU ijk(-’£}  ti(-)  ~  Eijk(-’^)  d? 
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in  which  the  stress  field  is  now  expressed  in  terms  of  boundary  tractions, 

boundary  displacements,  and  body  forces.  Since  the  compatibility  relation 

(Equation  (2))  involves  derivatives  of  displacement,  the  kernel 

2 

functions  H  and  E  have  singularities  of  orders  1/r  and  1/r  , 
respective ly . 

The  DEEM  is  formulated  by  the  numerical  approximation  of  Equations 
(10).  Equation  (10a)  is  used  to  obtain  the  unknown  boundary  values,  and 
then  both  equations  allow  the  calculation  of  internal  responses.  The 
two  approximations  made  to  obtain  the  unknown  boundary  values  are:  (1) 
integrating  in  a  piecewise  manner,  and  (2)  solving  the  equation  in  a 
boundary  collocation  sense.  That  is,  the  integration  is  subdivided  over 
boundary  elements  and  domain  cells;  and  the  integral  equations  are  applied 
to  a  discrete  number  of  points  on  the  boundary. 

The  early  applications  of  the  method  used  analytical  integrations 
over  elements  which  could  model  constant  or  linear  distributions  of 
boundary  and  domain  values.  More  recent  applications  apply  isoparametric 
element  concepts  common  to  the  finite  element  method  (Ref  32).  However, 
unlike  the  FEM  all  quantities  are  interpolated  at  the  same  order  (i.e., 
linear  boundary  elements  interpolate  geometry,  displacements,  and 
tractions  linearly).  The  boundary  distribution  of  any  field  variable  J 
on  a  single  element  is  written  as 


(.(n)  =  l 


(ri)  y, 


where  a  is  the  node  index;  N„  is  the  number  of  element  nodes;  ri  is  the 

LN 

normalized  local  curvalinear  coordinate;  If,  are  nodal  values  of  the 

ia 

field  variable;  and  N^ln)  are  the  appropriate  polynomial  shape  functions 
for  the  element.  A  typical  quadratic  element  (N^  =  3)  is  shown  in 
Figure  4.  The  shape  functions  for  this  element  are  given  by 


N i ( n )  =  |  (n2  -  n) 


N2(n)  =  (l  -  n  ) 
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N3(n)  =  \  (n2  +  n) 
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The  kernal  functions  in  the  integral  equations  are  singular  and  require 
special  treatment  when  considering  the  response  at  points  on  (Ref  1)  or 
near  (Ref  30)  the  boundary.  The  corresponding  numerical  details  are 
omitted  here  for  brevity.  While  the  use  of  the  fundamental  solution 
ensures  the  satisfaction  of  the  governing  differential  equations,  the 
boundary  collocation  ensures  violation  of  the  boundary  conditions.  This 
point  is  discussed  in  a  later  section  with  respect  to  a  coupled  solution 
approach . 

The  piecewise  integration  and  boundary  collocation  allow  Equation 
(10a)  to  be  approximated  by  a  system  of  linear  equations  in  terms  of  the 
nodal  boundary  values  and  body  forces  as 

Gt-Fu  +  G'<j>  =  0  (13) 

(for  detail  see  Reference  1  or  2).  The  body  forces  are  specified  and 
thus  the  product  with  G1  gives  a  known  vector.  For  a  mixed  boundary 
value  problem  the  vectors  of  nodal  tractions  and  displacements  contain 
known  and  unknown  values.  To  solve  for  the  unknown  nodal  values  we 
partition  the  system  as 
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+  G'  =  0 


and  then  collect  the  unknown  nodal  values  to  give 


Si  -El 


-G.  F 
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The  unknown  nodal  values  are  obtained  by  solving  the  above  system.  Then 
Somigliana's  identities  can  be  used  to  obtain  the  desired  internal 
responses . 


Overview  of  Indirect  Boundary  Element  Method 
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The  indirect  boundary  element  method  is  probably  the  second  most 
common  integral  equation  method.  Like  the  direct  method  it  also  has  its 
origins  in  classical  potential  theory  (Ref  35).  Single-  and  double¬ 
layer  potentials  were  used  in  the  theory  of  classical  electrodynamics  to 
express  boundary  value  problems  as  integral  equations.  The  most  advanced 
implementation  of  the  method  appears  to  remain  in  its  application  to 
electromagnetic  field  problems  (Refs  6,  17,  and  36).  The  indirect 
method  appears  to  have  first  been  applied  to  elastostatics  by  Massonnet 
(Ref  37).  The  numerical  development  of  IBEM,  for  elasticity  problems, 
paralleled  that  of  the  DBEM  until  the  mid  1970s.  In  recent  years  it  has 
not  been  developed  as  extensively  as  the  direct  method;  however,  its 
physically  meaningful  formulation  provides  insight  to  both  methods 
(Ref  31). 

The  integral  equations  on  which  the  IBEM  is  based  are  a  single¬ 
layer  potential  statement  of  the  boundary  value  problem.  The  domain 
ft  is  embedded  in  an  infinite  plane  as  shown  in  Figure  5.  For  elasticity 
applications  the  single- layer  source  corresponds  to  a  vector  of 
artificial  tractions,  P^(£)  (£  £  O.  In  this  formulation  we  seek  the 
boundary  distribution  of  artificial  tractions  which  satisfy  the 
prescribed  boundary  conditions.  These  tractions  are  artificial  in  the 
sense  that  they  only  exist  because  the  domain  ft  has  been  embedded  in  an 
infinite  plane.  They  represent  an  intermediate  step  in  the  formulation 
and  do  not  correspond  to  the  actual  boundary  tractions.  Artificial  trac¬ 
tions  and  body  forces  can  be  expressed  as  a  "continuous"  distribution  of 
e^(£) .  Thus  the  displacement,  stress  and  traction  fields  are  given  by 
integrating  Equations  (7a),  (8a),  and  (9a)  respectively  as 
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ti(x)  =  I  Fik(x,£)  Pk(|)  d%  +  f  Fik(x,z)  <^k(z)  dz 
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Since  the  solution  is  expressed  as  a  superposition  of  the  fundamental 
solution,  the  governing  differential  equations  are  satisfied  over  the 
entire  plane  including  the  domain  S2. 

To  determine  the  distribution  of  the  artificial  tractions  we  bring 
the  field  point  x  to  the  boundary  T  and  enforce  the  boundary  conditions, 
Equation  (6).  The  resulting  integral  expressions  are  given  as 
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u.(x)  =  I  G.k(x,p  Pk(|)  d£  +|  Gik(x,z)  4>k(z)  dx  x  e  (19) 


t.(x)  =  j  FiR(x,p  Pk(£)  d|  +  J  Flk(x,z)  *k(z)  dz 


x  e  rT  (20) 
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Equation  (19)  is  regular  upon  integration  while  Equation  (.20)  must  be 
interpreted  in  a  Cauchy  principal  value  sense  and  is  thus  written  as 
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A  tangent  line  is  assumed  through  x  and  the  sign  on  the  first  term 
depends  on  the  orientation  of  the  element  with  respect  to  fi. 

The  IBF.M  is  formulated  by  the  numerical  approximation  of  Equations 
(16)  through  (21).  Equations  (19)  and  (21)  are  used  to  determine  the 
unknown  artificial  tractions,  and  then  Equations  (16)  through  (18)  allow 
the  calculation  of  internal  responses.  As  with  the  DBEM,  the  two 
approximations  made  to  obtain  the  unknowns  are:  (1)  integrating  in  a 
piecewise  manner  and  (2)  solving  the  integral  equations  in  a  weighted 
residual  sense  (on  the  boundary).  The  integration  is  subdivided  over 
boundary  elements  and  domain  cells.  Normally  the  integral  equations  are 
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satisfied  in  a  collocation  sense;  however,  Lean  et.al  (Ref  36)  report 
improved  accuracy  by  a  Galerkin  approximation.  This  study  is  limited  to 
a  collocation  approximation  of  the  equations.  Equations  (19)  and  (21) 


are  then  approximated  as 


u  =  G  P  +  G' 


t  =  F  P  +  F' 


(22a) 


(22b) 


where  u  and  t  are  values  of  boundary  displacements  and  tractions  at  col¬ 
location  points.  For  continuous  elements  the  collocation  points  correspond 
to  element  node  points.  The  coefficient  matrices  G,  G' ,  F  and  F'  are 
obtained  by  integrations  of  Equations  (19)  and  (21)  with  respect  to  the 
appropriate  shape  functions.  P  is  the  vector  of  unknown  nodal  artificial 
traction  values,  and  the  last  term  provides  the  effect  of  the  body  forces 
at  each  collocation  point.  For  details  see  References  1,  2,  31,  and  38. 

For  a  mixed  boundary  value  problem  Equations  (22a)  and  (22b)  are 
applied  to  and  f^,  respectively,  giving 


The  nodal  artificial  tractions  are  obtained  by  solving  the  above  system. 

Then  Equations  (16)  through  (18)  can  be  used  to  obtain  the  desired  internal 

responses.  As  when  developing  the  system  of  equations,  boundary  values 

are  obtained  by  letting  x  go  to  the  boundary  (as  Equations  (19)  and  (21)). 

Consider  a  few  of  the  differences  between  the  two  BE  formulations. 

Let  N(,p  be  the  number  of  boundary  collocation  points.  Ignoring  the 

integrations  associated  with  body  forces,  the  DBEM  (Equation  15)  and 

2  2 

IBEM  (Equation  23)  require  the  calculation  of  2(2N^p)  and  (2N^,p) 

coefficients,  respectively.  However,  in  solving  the  system  of 

equations,  the  DBEM  yields  the  unknown  boundary  values  directly  while 

the  IBF.M  yields  the  artificial  traction  distribution.  To  calculate  all 

2 

the  unknown  boundary  values  the  IBEM  requires  additional  (2N„D) 


coefficient  calculations.  Though  the  total  number  of  coefficients,  and 
thus  integrations,  required  to  obtain  the  unknown  boundary  values  are 
equal,  the  computational  effort  is  not.  The  coefficient  matrices  of  the 
DBF.M  are  calculated  "simultaneously"  and  thus  the  "overhead  calculations" 
(e.g.,  calculation  of  the  jacobian)  associated  with  the  numerical  integra¬ 
tions  are  done  once.  Since  these  "overhead  calculations"  are  performed 
twice  by  the  IBEM,  it  requires  more  effort.  However,  this  comparison 
has  assumed  the  analyst  needs  all  the  boundary  values.  They  are  required 
by  the  DBF.M  to  calculate  internal  responses  by  the  Somigliana  identities; 
they  are  not  needed  by  the  IBEM  for  internal  response  calculation  and 
thus  the  analyst  can  be  selective. 

The  calculation  of  internal  responses  also  differs  between  the 
methods.  The  DBEM  obtains  internal  responses  by  the  Somigliana  identi¬ 
ties,  Equation  (10),  which  integrate  the  effects  of  both  the  boundary 
tractions  and  displacements.  The  IBEM  only  integrates  the  effects  of 
artificial  tractions.  Equations  (16)  through  (18).  As  in  developing  the 
system  of  equations,  though  the  effort  is  not  doubled  for  the  DBEM,  it 
is  considerably  more.  Another  important  difference  in  the  internal 
response  calculations  is  the  order  of  singularity  of  the  kernal 
functions.  For  two  dimensional  problems  the  methods  have  the  following 
orders  of  singularity: 


Displacement 

o( In  r )  and  o( 1/r) 
o( In  r) 


Strain,  Stress,  Traction 

o( 1/r)  and  o( 1/r^) 
o( 1/r) 


Thus  the  DBEM  must  deal  with  stronger  singularities  and  correspondingly 
more  difficult  integrations.  This  problem  is  most  severe  in  the  near 
boundary  region  (Refs  30  and  31).  In  fairness,  the  indirect  method  is 
not  without  its  problems  in  calculating  internal  responses.  Though  the 
integration  effort  is  reduced  in  the  IBEM,  often  the  accuracy  is  too.  A 
problem  that  is  inherent  to  the  IBEM  is  associated  with  loading  and 
geometric  discontinuities.  In  these  areas,  the  artificial  tractions 
experience  very  high  gradients  even  though  boundary  conditions  may  be 
very  well  behaved.  Thus  unless  these  areas  receive  special  treatment  in 
the  numerical  formulation  or  modeling,  accuracy  is  locally  very  poor. 
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For  coupled  solution  approaches,  particularly  for  infinite  domain 
problems,  it  is  not  easy  to  determine  which  BEM  would  be  the  most 
effective.  Both  methods  must  calculate  two  coefficient  matrices  which 
involve  the  same  orders  of  singularity.  In  addition,  discontinuities 
can  often  be  avoided,  eliminating  them  as  a  factor.  In  the  following 
section  the  IBEM  is  used  to  illustrate  a  direct  derivation  of  a 
stiffness  matrix.  The  ideas  can  be  applied  to  the  DBEM  but  are 
explained  with  the  IBEM  due  to  its  conceptual  simplicity. 

A  Direct  Derivation  of  an  IBEM  Stiffness  Matrix 

This  section  presents  a  physically  intuitive  derivation  of  an  IBEM 
stiffness  matrix.  The  relation  between  nodal  displacements  and  tractions 
is  identical  to  the  relation  obtained  directly  by  Kelly  et  al.  (Ref  11). 
The  physical  approach  provides  additional  insight  to  the  stiffness 
formulat ion . 

A  stiffness  matrix  relates  nodal  displacements  to  generalized  nodal 
forces  in  the  form 


K  u  =  f  (24) 

where  K,  u  and  f  are  the  stiffness  matrix,  displacement  vector,  and 

generalized  force  vector,  respectively.  This  basic  definition  and  a 

numerical  solution  provided  by  the  IBEM  allow  a  stiffness  matrix  to  be 

calculated  for  a  BEM  region.  Consider  the  displacement  boundary  value 

problem  shown  in  Figure  6.  The  boundary,  T,  is  subdivided  into  n„ 

h 

boundary  elements.  The  elements  are  shown  as  isoparametric  quadratic 
elements  for  i 1  lust  rat  ion .  All  displacements  are  zero  except  the 
d  o  f  which  has  a  unit  displacement.  We  can  now  use  the  IBEM  to  solve 
the  displacement  boundary  value  problem.  As  given  by  Equation  (22a) 
with  body  forces  omitted,  the  unknown  nodal  artificial  boundary  trac¬ 
tions  are  related  to  the  known  nodal  displacements  by 


u  =  G  P 
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This  system  of  equations  can  be  solved  for  the  nodal  artificial  traction 
values  which  can  then  be  used  to  obtain  the  unknown  real  boundary  trac¬ 
tions  by  an  approximation  of  Equation  (18).  The  nodal  tractions  can  be 
expressed  in  matrix  form  as 

t  =  F  P  (26) 

The  distribution  of  the  traction  along  the  boundary  is  then  approximated 
in  terms  of  the  locally  based  shape  functions*  as 

t(x)  =  M(x)  t  (27) 

which  by  incorporating  Equation  (26)  becomes 

f ( x )  =  M( x )  F  P  (28) 

With  the  distribution  of  boundary  tractions,  generalized  nodal  forces 
(on  the  boundary)  can  now  be  obtained.  The  nodal  forces  are  determined 
by  "measuring"  the  work  done  by  the  tractions  due  to  a  series  of  virtual 
displacements  -  displacement  shape  functions*.  The  generalized  forces 
are  given  by 

f  =  f  NT(x)  t(x)  dx  (29) 

r 

For  the  prescribed  displacements  u,  the  generalized  forces  (Equation 
(29))  are  the  column  of  a  stiffness  matrix  for  the  domain  fl.  To 

completely  calculate  the  stiffness  matrix  we  must  individually  subject 
each  remaining  nodal  degree  of  freedom  to  a  unit  perturbation,  solve  the 
corresponding  displacement  boundary  value  problem,  and  then  determine 
the  generalized  forces  by  Equation  (29)  to  obtain  the  corresponding 
column  of  K. 

We  now  seek  to  express  the  above  procedure  in  equation  form.  Using 
indicia]  notation,  the  nodal  displacement  vector  for  the  J ^ 
perturbation  is  given  by 

*In  this  context  the  shape  functions  are  locally  based  but  defined 
globally.  That  is,  they  are  nonzero  on  one  or  two  elements  but  defined 
ori  all  f . 


-J 
u  . 
1 


1,  i  =  J 
0.  i  #  J 


(30) 


where  i  ranges  from  1  to  the  number  of  degrees  of  freedom  (n).  Each  u 
is  the  coordinate  vector  in  the  n  dimensional  space.  Combining  all 
of  the  vectors  into  a  single  matrix  (i.e.,  the  superscript  J  becomes  an 
index  j)  allows  Equation  (25)  to  be  rewritten  as 


u  .  . 
'  1 


ik 


( 31) 


By  Equation  (30),  u  is  simply  the  identity  matrix.  Thus  we  can  easily 
solve  for  t  lie  artificial  tractions  as 


P.  . 

1.1 


(32) 


This  formulation  provides  insight  to  the  character  of  G  * ;  each  column 
of  G  '  is  the  vector  of  artificial  tractions  resulting  from  a  unit 
perturbation  of  the  corresponding  d  o  f.  G  is  of  full  column  rank  and 
thus  invertable  if  the  collocation  point  positions  are  unique. 

Inserting  this  result  into  Equations  (28)  and  (29)  yields  the  stiffness 
matrix  as 


K.  .  -  /  N.  .  ( x )  M.  ,  ( x )  dx  F.  G  . 
ij  J  ki  -  kl  -  -  lm  mj 

r 

or  in  matrix  notation 


(33) 


K  =  J  N^(x)  Mix)  dx  F  G_1  (34) 

r 

Kelly  et.  al.  (Ref  11)  obtained  the  same  form  of  solution  for  the  potential 
problem  by  eliminating  the  source  density.  For  the  e 1  as t os  tat l cs  problem 
this  is  equivalent  to  eliminating  P  from  Equations  (25)  and  (26).  With 
their  approach  we  see.  that  the  product  of  F  and  G  '  provides  the  relation 
between  nodal  tractions  and  displacements.  We  can  write  Equation  (34) 
in  an  abbreviated  form  as 


2! 


. .  .  . 


K  =  C  E 


where 


C  =  j  N^(x)  M(x)  dx 


E  =  F  G  *  (37) 

Brebbia's  direct  approach  (Ref  21)  to  obtaining  a  stiffness  matrix  for 
the  DBEM  also  has  the  same  form.  For  the  DBEM  however,  E  is  given  by 

E  =  g"1  F  (38) 

wliore  F  Is  not  the  same  as  obtained  for  the  IBEM.  The  DBEM  uses  the 
same  fundamental  solution;  however  the  roles  of  the  source  and  field 
points  and  the  roles  of  the  indices  are  reversed  in  obtaining  the 
Somigliana  identities. 

As  previously  mentioned,  direct  formulations  of  stiffness  matrices 
do  not  initially  yield  a  symmetric  matrix.  In  addition  both  the  direct 
and  variational  formulations  produce  stiffness  matrices  which  violate 
equilibrium.  Hartman  (Ref  23)  attributed  these  problems  to  the  approxi¬ 
mations  made  in  obtaining  the  traction-displacement  relationship,  Equa¬ 
tions  (37)  and  (38),  for  the  IBEM  and  DBEM,  respectively.  For  the  DBEM 
the  equilibrium  problem  can  be  addressed  when  formulating  E,  liquation 
(13)  is  augmented  with  two  additional  equations  and  corresponding  Lagrange 
multipliers,  which  force  the  traction  distribution  to  "identically" 
satisfy  equilibrium  (Ret  11).  Most  approaches  to  these  problems  modify 
the  stiffness  matrix  rather  than  address  the  source  of  the  problem,  the 
tract  ion -displacement  relationship  (Ref  20). 

Discontinuity  of  Boundary  Element  Methods 

Regardless  of  the  formulation  there  is  a  characteristic  of  BE 
formulated  stiffness  matrices  and  BE  solutions  in  general  which  seems  to 
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be  ignored  in  the  literature;  BE  formulations  are  inherently 
incompatible  at  the  interface  of  homogeneous  regions. 

Derivations  of  BE  formulations  usually  indicate  that  the  governing 
differential  equations  are  satisfied  exactly  in  the  domain  while  the 
boundary  conditions  are  only  approximately  satisfied.  Often  this  is 
interpreted  as  meaning  the  boundary  element  shape  functions  constrain 
the  displacement  field.  Actually  the  meaning  has  greater  depth.  This 
interpretation  is  often  a  consequence  of  invalid  analogies  between  BEMs 
and  FEMs .  The  methods  are  theoretically  related  and  certainly  the  BEM 
has  borrowed  much  technology  from  the  FEM;  however,  the  accepted  use  of 
the  name  "boundary  element  methods"  for  boundary  integral  equation  tech¬ 
niques  is  misleading.  This  interpretation  is  usually  not  enlightened  by 
numerical  experience  since  many  BE  programs  can  not  accurately  calculate 
responses  in  the  near  boundary  region. 

The  shape  functions  used  for  the  so  called  "boundary  elements"  do 
not  constrain  the  displacement  field  but  merely  offer  an  approximation 
of  the  boundary  values  for  integration.  If  they  did  constrain  the 
displacement  field,  a  single  BE  region  comprised  of  four  linear  elements 
could  pass  the  FE  "patch  test"  (Ref  39);  such  is  not  the  case.  Typically 
the  integral  equations  are  satisfied  in  a  collocation  sense  (i.e.,  at 
discrete  points  along  the  boundary)  and  thus  the  boundary  conditions  are 
satisfied  at  the  collocation  points  (which  coincide  with  nodal  points 
for  continuous  element  formulations) .  In  between  these  points,  satis¬ 
faction  of  boundary  conditions  can  be  greatly  in  error.  This  is  graph¬ 
ically  illustrated  in  our  previous  report  (Ref  31). 

Similar  to  the  displacement-based  FEM's  satisfaction  of  equilibrium 
equations  in  a  nodal  sense,  the  BEM  only  satisfies  boundary  conditions 
in  a  nodal  sense.  Many  researchers  indicate  that  compatibility  between 
the  two  methods  is  easily  satisfied  by  using  the  same  shape  functions 
(Refs  1,  13,  21,  22,  and  24);  I  disagree.  However,  discontinuity  does 
not  imply  poor  results.  On  the  contrary  many  incompatible  finite 
elements  exhibit  improved  performance,  and  many  accurate  BE  solutions 
have  been  reported.  Recognition  of  this  basic  incompatible  behavior 
could  be  of  practical  use.  It  affects  the  convergence  properties  of  the 
method  and  can  potentially  explain  other  solution  characteristics.  One 


very  useful  application  could  be  in  mesh  refinement.  Whether  the 
refinement  is  a  manual  or  an  adaptive  process,  the  variation  in  computed 
responses  between  collocation  points  can  provide  a  measure  on  which  to 
base  the  refinement. 


DISCUSSION  OF  NUMERICAL  IMPLEMENTATION 

This  section  discusses  the  numerical  implementation  of  the  IBEM 
stiffness  matrix  deve^ped  above.  Algorithms,  in  the  form  of  pseudo¬ 
code  outlines,  provide  detail  on  key  aspects  of  the  implementation. 

These  algorithms  have  been  implemented  in  a  research  code.  The 
numerical  studies  are  incomplete  but  will  be  presented  in  a  follow-up 
report.  The  scope  of  this  discussion  is  limited  to  the  calculation  of 
the  stiffness  matrix;  assembly  of  element  stiffness  matrices  is  well 
documented  (Ref  39).  Both  the  IBEM  and  DBEM  are  addressed  below  since 
their  stiffness  matrices  Equation  (35)  can  have  the  same  form.  Rudolphi 
(Ref  22)  provides  a  general  outline  for  the  calculation  of  a  "stiffness" 
matrix  including  the  BE  coefficient  matrix  calculations .  His  paper 
deals  principally  with  the  potential  problem  using  the  DBEM.  Though 
Rudolphi's  work  was  not  referenced  for  the  implementation  aspects  of 
this  work,  the  main  steps  in  the  calculations  are  independent  of  the 
application  and  BE  formulation. 

Due  to  the  complexity  of  BE  and  FE  software  systems  emphasis  must 
be  placed  on  the  modular  design  of  the  coupled  software  system.  In  the 
following  discussion  I  assume  BE  and  FE  systems  exist  and  require  modi¬ 
fication.  The  stiffness  matrix  calculation  can  be  organized  as  three 
tasks : 

(1)  Calculation  of  the  coefficient  matrices  (G  and  F)  and  C 

(2)  Preliminary  K  calculation,  inverse  and  matrix  product 
calculations  as  shown  in  Equation  (34) 


(3)  Optional  adjustments  in  K  to  satisfy  equilibrium  and 
symmetry 


The  first  task  requires  modification  of  the  BE  system.  The  second 
and  third  tasks  can  be  added  to  the  BE  system,  comprise  a  separate 
module,  or  be  incorporated  into  an  element  routine  of  the  FE  system.  In 
my  approach  I  prefer  to  "weakly  couple"  the  software  systems  so  that 
they  can  still  serve  as  individual  research  and  analysis  tools.  By  this 
approach  the  first  task  is  completely  performed  by  the  modified  BE 
system.  The  second  and  third  tasks  are  performed  by  a  separate  program 
module  which  except  for  the  three  component  matrices  of  K  (i.e., 

C,  F,  G,)  requires  no  access  to  the  BE  data.  An  element  routine  is  then 
added  to  the  FE  program  which  reads  data  defining  the  "super-element" 
connectivity  and  assembles  the  BF,  K  into  the  global  system  at  the  time 
of  assembly.  Solution  of  the  FE  system  provides  nodal  displacements 
from  which  stresses  can  be  calculated  in  both  the  FE  and  BE  domains. 

When  boundary  tractions  or  internal  responses  are  required  in  the 
BE  region,  the  displacements  associated  with  BE  nodes  must  be  retained. 

For  either  BE  formulation,  boundary  tractions  can  be  obtained  by  the 
equation:  t  =  E  u.  For  the  DBEM  the  known  boundary  values  in  combination 

with  Somigliana's  identities  (Equation  (10))  provide  internal  responses. 

For  the  I BEM  the  tractions  or  displacements  could  be.  used  to  solve  a 
boundary  value  problem  for  the  artificial  tractions.  Alternatively  G  * 
could  be  saved  during  the  K  calculations  giving  P  =  G  1  u.  The  effect 
of  the  artificial  tractions  is  then  integrated  to  obtain  internal  responses 
according  to  Equations  (16)  through  (IB).  The  remainder  of  this 
discussion  concentrates  on  the  calculation  of  the  stiffness  matrix. 

For  smaller  problems  the  calculation  of  the  matrices  comprising  K 
dominates  the  numerical  effort.  The  coefficient  matrices  are  inherently 
full  arid  their  formulation  requires  extensive  numerical  integration. 

The  cost  of  their  calculation  would  appear  to  increase  with  the  square 
of  the  number  of  boundary  elements.  There  is  an  element  "overhead  cost" 
which  follows  this  trend;  however  for  a  program  which  uses  a  variable 
orrler  integration  scheme  the  cost  of  the  average  element  integration 
decreases.  This  is  due  to  the  reduction  in  the  order  of  integration 
with  a  decrease  in  the  ratio  of  e.lement  length  to  collocation  point 
distance  (Ref  31  and  32). 
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Calculation  of  E  requires  an  IBEM  program  to  generate  both  F  and  G 
for  a  single  homogeneous  domain.  These  are  the  same  equations  required 
to  combine  homogeneous  BE  regions  by  approximate  satisfaction  of  equi¬ 
librium  and  continuity  along  their  interface.  Thus  many  existing  BE 
programs  already  have  this  capability.  The  calculation  of  two  coeffi¬ 
cient  matrices  does  not  double  the  numerical  effort  because  much  of  the 
overhead  in  the  numerical  integrations  is  common  to  both  types  of  co¬ 
efficients  . 

In  the  DBEM  both  F  and  G  are  always  calculated.  The  only  modifica¬ 
tion  which  might  be  required,  depending  on  the  program  design,  is  the 
retention  of  all  coefficients.  Some  implementations  immediately  multiply 
known  boundary  values  by  corresponding  coefficients  and  sum  the  terms 
into  the  known  vector  to  form  the  system  of  equations  (see  Equation  15). 

Calculating  the  integration  of  the  shape  function  product  matrix, 

C,  according  to  Equation  (36),  is  effectively  performed  at  the  element 
level.  Because  of  the  locally  based  nature  of  the  shape  functions,  C  is 
very  sparse.  For  quadratic  elements,  rows  associated  with  a  mid-node 
have  three  nonzero  terms,  and  those  associated  with  an  extreme  node  have 
five  nonzero  terms.  In  this  implementation  tractions  and  displacements 
are  interpolated  at  the  same  order  (quadratic),  and  there  are  no 
provisions  for  traction  discontinuities.  As  a  result  M(x)  =  N(x)  and 
therefore  C  is  symmetric.  Not  providing  for  traction  discontinuities 
does  place  geometrical  constraints  on  the  BE  region.  However  for  BE 
regions  comprised  of  several  elements  and  few  corners,  Rudolphi  (Ref  22) 
concludes  that  traction  continuity  has  a  negligible  effect  on  the  accuracy 
For  appl ications  where  the  BEM  region  is  used  as  an  "infinite  element" 
the  analyst  determines  the  shape  of  the  FEM/BEM  interface  and  thus  geo¬ 
metric  discontinuities  pose  no  problem. 

For  a  vector  boundary  value  problem  such  as  the  elasticity  problem, 
the  two  rows  of  C  corresponding  to  a  common  node  have  identical  terms 
but  differ  in  their  position,  column.  Except  for  the  diagonal  terms  in 
C  each  term  calculated  at  the  element  level  is  complete  (i.e.,  no  sum¬ 
mation  is  required)  and  could  be  formally  assembled  into  four  positions 
of  C.  Only  the  diagonal  terms  of  C  which  correspond  to  shape  functions 
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spanning  two  elements  require  the  addition  of  a  second  term.  Thus  it  is 
very  efficient  to  both  calculate  C  and  perform  its  multiplication  with  E 
at  an  element  level.  There  is  a  maximum  of  six  unique  terms  for  a  single 
element.  Execution  and  storage  requirements  associated  with  C  are  trivial 
compared  to  F  and  G  which  are  full  coefficient  matrices. 

The  algorithm  below  outlines  the  calculation  of  C  at  the  element 
level.  The  research  code  is  written  in  the  Modula-2  language  (Ref  40), 
which  is  a  descendent  of  Pascal  and  Modula.  Those  familiar  with  Pascal 
will  recognize  the  dialect  below.  I  have  attempted  to  remove  enough 
syntactical  idiosyncrasies  from  the  code  to  permit  communication.  Calls 
to  Modula-2  procedures  (analogous  to  FORTRAN  subroutines)  are  supplemented 
with  short  descriptions  of  their  purpose.  Comments  which  do  not  replace 
actual  code  are  enclosed  in  (*  *)'s.  General  file  operations  and  variable 
definitions  are  omitted.  Descriptive  variable  names  supplemented  by 
comments  are  used  to  define  variables.  Many  modern  languages  provide 
for  user  defined  data  types.  Nnode  and  Mnode  are  variables  which  can 
have  the  enumerated  values  of  (a,c,b),  the  node  indices  listed  by  extreme 
nodes  first.  This  corresponds  to  a  of  Figure  4  taking  the  values  of 
(1,3,2)  respectively.  These  variables,  Nnode  and  Mnode,  serve  as  nodal 
indices  to  arrays  and  promote  internal  documentation.  More  complicated 
data  types  such  as  records  have  been  converted  to  arrays  to  simplify  the 


PROCEDURE  Genera teShapeProductMatr ix 

(ContNode,  (*  array  of  continuous  element  node  numbers 

indexed  as  (element_number ,node_index)  *) 
numeltotal,  (*  total  number  of  continuous  elements  *) 
ShapelntOrder )  (*  numerical  integration  order  *) 

BEGIN  (*  the  GenerateShapeProductMatrix  procedure  *) 

FOR  elnum  =  1  TO  numeltotal  DO  (*  for  each  element  *) 

(*  determine  shape  function  (N)  &  weight-Jacobi an  product 
values  at  each  integration  point  *) 

FOR  Nnode  =  a  TO  b  DO  (*  an  index  on  the  node  *) 

Write  the  equation  number  for  node  ContNode[ elnum, Nnode]  to 
the  NMfile 

END  (*  the  Nnode  loop  *) 

FOR  intpt  =  1  TO  ShapelntOrder  DO  (*  each  integration  point  *) 
position  =  GAUSSpt[ ShapelntOrder, intpt]  (*  obtain  the 
position  of  the  integration  point  *) 
SIlAPEfunctions(position,Nvalues, intpt)  (*  calculate  the  shape 
function  values  at  "position"  and  save  in  Nvalues 
indexed  as  [ intpt , node_index]  *) 

Calculate  the  Jacobian  at  "position" 

WGTdetJ] intpt]  =  GAUSSwgt[ ShapelntOrder, intpt]*J  (*  calculate 
the  integration  weight  Jacobian  product  at  each  integration 
point  *) 

END  (*  the  integration  point  loop  *) 

FOR  Nnode  =  a  TO  b  DO  (*  for  each  nonzero  shape  function  *) 

FOR  Mnode  =  Nnode  TO  b  DO  (*  for  each  unique  combination  of 
shape  function,  i.e.,  N=M  thus  symmetry  is  considered  *) 
integral  =0.0  (*  initialize  integration  *) 

FOR  intpt  =  1  TO  ShapelntOrder  DO  (*  each  int.  pt  *) 
integral  =  integral  +  Nva lues [ intpt, Nnode]* 

Nvaluesf intpt, Mnode ]*WGTdetJ[ intpt] 

END  (*  the  integration  point  loop  *) 

Write  the  integral  value  to  the  NMfile 
END  (*  the  Mnode,  M  shape  function,  loop  *) 

END  (,*  the  Nnode,  N  shape  function,  loop  *) 

END  (*  the  element  loop  *) 

END  (*  the  GenerateShapeProductMatrix  procedure  *) 
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The  six  integration  results  are  written  to  the  NMfile  in  the  order: 
a-a,  a-c,  a-b,  c-c,  c-b,  and  b-b,  where  the  letters  indicate  the  node 
index  of  each  shape  function  in  the  product. 

With  all  of  the  component  matrices  for  K  calculated,  the  second 
task  of  calculating  an  initial  nonsymmetrical  stiffness  matrix  can  pro¬ 
ceed.  The  third  task  of  satisfying  equilibrium  and  obtaining  a  symmetric 
form  can  be  combined  with  the  matrix  operations  of  the  second  task.  The 
details  of  integrating  the  two  tasks  are  highly  dependent  on  the  method 
used  to  obtain  symmetry  and  equilibrium.  For  this  research  code,  the 
tasks  are  independent,  so  different  methods  can  be  investigated.  This 
also  allows  the  symmetry  of  the  original  stiffness  matrix  to  be  examined. 

The  initial  steps  required  in  the  calculation  of  K  differ  in  the 
two  BE  formulations.  This  difference  results  from  the  reversal  of  F  and 

G  ^  in  the  traction-displacement  relations,  Equations  (37)  and  (38).  As 

3 

a  result  of  this  difference  the  IBEM  is  burdened  with  n  additional  multi- 

2 

plications  and  (n-l)n  additions  to  obtain  E.  Additionally  the  DBEM  has 
reduced  memory  requirements.  F  can  be  read  from  disk  one  column  at  a 
time.  Each  column  can  immediately  be  used  to  calculate  the  correspond¬ 
ing  column  of  K. 

(V 

As  the  order  of  the  system  of  equations  increases,  the  calculation 

of  the  inverse  in  Equation  (37)  or  (38)  becomes  the  most  costly  step. 

3 

For  Gauss  elimination  the  cost  increases  as  n  .  A  common  approach  (Ref 
22  and  24)  to  reducing  the  cost  of  calculating  the  inverse  is  to  sub¬ 
divide  the  single  BE  region  into  several  BE  regions,  providing  many  small 
stiffness  matrices  instead  of  one  large  stiffness  matrix.  This  subdivi¬ 
sion  introduces  additional  boundary  elements,  and  thus  more  equations, 
but  the  equations  corresponding  to  the  complete  BE  domain  are  now  block 
banded.  Mitsui  et  al.  (Ref  24)  suggests  that  with  regard  to  efficiency 
there  is  an  optimal  degree  of  subdivision.  Subdividing  the  BE  region 
may  also  require  that  traction  discontinuities  and  the  violation  of  equi¬ 
librium  be  treated  more  rigorously  since  the  occurrence  of  geometric 
discontinuities  and  a  relatively  small  number  of  elements  would  be  more 
1  i  ke  1  y  . 
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For  both  BE  formulations,  a  column  of  E  is  immediately  used  to 
calculate  the  corresponding  column  of  K.  The  above  algorithm  calculates 
C  at  the  element  level  and  writes  the  values  to  NMfile.  An  algorithm  to 
calculate  an  initial  K,  given  the  component  matrices,  is  shown  below. 
Values  of  C  are  read  from  NMfile,  consistent  with  the  previous 
algorithm.  In  this  implementation  the  pseudo-code  given  below  is  a 
separate  program  (a  Modula-2  program  module). 


VV  ' - 


MODULE  Stif fnessMatrix 

BEGIN  (*  the  Stif fnessMatrix  module  *) 

ReadCmatrix(EqNum,C,numeltotal .order)  (*  read  C:  EqNum  contains 
the  equation  numbers  for  each  node  indexed  as 
[ element_number ,node_indexj  and  C  contains  the  integration  of 
the  shape  function  products  indexed  as  [ element_number , 0. . 5 ]  *) 
ReadGmatrix  (*  read  G  *) 

ReadFmatrix  (*  read  F:  For  the  Modula-2  implementation  these  two 
procedures  cause  G  and  F  to  be  read  by  modules  envolved  in  the 
calculation  of  E.  This  module  does  not  have  access  to  the  data 
of  G  and  F.  *) 

(*  Calculate  the  stiffness  matrix.  K  *) 

IV 

Initialize  K  to  zero 

FOR  col  =  1  TO  order  DO  (*  for  each  column  of  K  *) 

CalcEcol(col ,Ecol)  (*  calculates  a  single  column  of  E  stored  in 
Ecol  *) 

FOR  elnum  =  1  TO  numeltotal  DO  (*  for  each  element  *) 

Ccount  =  0  (*  initialize  the  index  for  the  C  array  *) 

FOR  nodel  =  a  TO  b  DO  (*  node  index  for  the  first  shape 
function  *) 

roweq  =  EqNum[ elnum, nodel ]  (*first  row  equation  for  nodel*) 
FOR  node2  =  nodel  TO  b  DO  (*  for  each  new  combination  *) 
coleq  =  EqNum [ elnum, node2]  (*  first  column  equation  for 
node2  *) 

Cval  =  C[ elnum, Ccount]  (*  extract  the  C  term  from  the 
element  array  *) 

FOR  eqincrement  =  0  TO  1  DO  (*  for  each  nodal  dof  *) 

Crow  =  roweq+eqincrement 
Ccol  =  coleq+eqincrement 

(*  calc,  contributions  in  both  symmetrical  terms  *) 
K[Crow,col]  =  K[Crow,col]  +  Ecol[Ccol]*Cval 
IF  roweq#coleq  THEN  (*  not  a  diagonal  term  *) 
K[Ccol,col]  =  K[Ccol,col]  +  Ecol[Crow]*Cval 
END 

END  (*  the  eqincrement  loop  *) 
increment  Ccount 
END  (*  the  node2  loop  *) 

END  (*  the  nodel  loop  *) 

END  (*  the  elnum  loop  *) 

END  (*  the  col  loop  *) 

Adjust  K  for  equilibrium  and  symmetry 

Write  the  upper  triangular  portion  of  K  to  disk 

END  (*  the  St  if fnessMatrix  module  *) 
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As  previously  noted,  due  to  symmetry  and  the  vector  nature  of  the 
elasticity  problem  each  element  C  value  (except  diagonal  terms)  has  four 
positions  in  the  assembled  C.  The  "eqincrement  loop"  above  multiplies 
each  element  C  value  by  four  terms  (two  for  diagonal  terms)  in  the  column 
of  E  and  assembles  the  products  into  the  corresponding  four  K  terms.  By 
this  approach  the  sparsity  of  C  is  exploited  and  the  matrix  is  never 
assembled. 

The  last  task  in  calculating  a  stiffness  matrix  is  adjusting  K  to 
satisfy  symmetry  and  equilibrium.  The  work  of  Tullberg  and  Bolteus 
(Ref  20)  provides  guidance,  based  on  numerical  studies,  for  both  of  these 
problems.  As  previously  mentioned,  numerical  studies  have  shown  the 
nonsymmetrical  K  obtained  by  the  direct  approach  to  be  more  accurate 
than  the  symmetrical  form  obtained  by  a  variational  approach.  However 
when  combining  the  matrix  with  an  existing  FE  system  a  symmetric  form  is 
usually  preferable.  An  exception  can  occur  if  the  FEM  is  used  to  model 
plasticity  governed  by  nonassociative  flow  rule;  in  this  case  the  result¬ 
ing  FF,  system  can  be.  nonsymmetr i c .  In  this  study  the  stiffness  matrix 
is  symmetrized  by  averaging  the  off-diagonal  terms. 

In  the  initial  implementation  of  the  research  code,  equilibrium 
considerations  have  not  been  included.  Our  principal  interest  is  in 
infinite  domain  problems  where  previously  mentioned  techniques  are  not 
applicable.  For  finite  domain  problems,  one  of  the  methods  implemented 
by  Tullberg  and  Bolteus  (Ref  20)  to  solve  the  equilibrium  problem  might 
be  improved  upon.  The.  equilibrium  error  for  a  given  direction  was  cor¬ 
rected  by  adding  Lite  average  error  to  each  term.  This  approach  did  not 
consider  the  relative  magnitude  of  each  term.  A  better  approximation 
might  be  obtained  by  basing  the  error  distribution  on  the  normalized 
magnitude  of  each  term  where  the  sum  of  absolute  values  of  the.  stiffness 
terms  is  the  normalizing  factor. 
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SUMMARY  ANT)  CONCLUSIONS 


Fundamental  differences  in  the  formulations  of  the  finite  element 
and  boundary  element  methods  result  in  corresponding  strengths  and 
weaknesses.  A  combination  of  the  methods  using  the  strengths  of  each 
may  allow  some  classes  of  problems  to  be  solved  more  effectively. 
Applications  with  infinite  or  semi-infinite  domains  and  fracture 
mechanics  problems  are  candidates  for  combined  solution  approaches. 

The  DBEM  and  IBEM  are  the  most  common  integral  equation  methods  in 
use.  They  both  have  their  origins  in  classical  potential  theory  and  the 
fundamental  singular  solution  is  a  key  component  of  their  derivations. 
For  both  methods  the  approximations  common  to  most  formulations  are: 

(1)  integrating  in  a  piecewise  manner,  and  (2)  solving  the  integral 
equations  in  a  collocation  sense.  The  collocation  approximation  limits 
the  satisfaction  of  boundary  conditions  to  discrete  points;  for  coupled 
solutions  interface  boundaries  will  be  inherently  incompatible 

Coupling  approaches  can  be  categorized  based  on  several  parameters. 
The  most  basic  classification  indicates  which  numerical  method  hosts  the 
coupling  (i.e.,  is  the  final  form  of  the  equations  "BE  like"  or  "FF. 
like").  This  study  addresses  FE-hosted  or  stiffness  couplings  which  can 
be  further  categorized  by  the  derivation  of  the  stiffness  matrix. 

Either  the  stiffness  matrix  is  obtained  directly  or  via  a  variational 
statement  of  the  problem.  In  both  cases  the  key  matrices  are  the  fully 
populated  coefficient  matrices  generated  by  the  BF.Ms  and  a  sparse  matrix 
obtained  by  the  integration  of  shape  function  products.  The  BF.M 
coefficient  matrices  are  manipulated  to  establish  a  relationship  between 
the  boundary  tractions  and  displacements;  this  relation  is  then  used 
directly  to  obtain  a  stiffness  matrix  or  substituted  into  a  boundary 
variational  statement. 

The  two  major  problems  with  the  BE  stiffness  matrix  are  its  lack  of 
symmetry  and  violation  of  equilibrium.  A  variational  formulation 
inherently  produces  a  symmetric  system  which  is  equivalent  to  taking  the 
symmetric  component  of  the  corresponding  d i red  1 y- formu 1 ated  stiffness 


n 
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matrix.  Investigators  have  offered  several  alternatives  for  dealing  with 
the  equilibrium  problem  ranging  from  the  introduction  of  Lagrange 
multipliers  which  enforce  equilibrium  to  ad  hoc  adjustments  of  the 
stiffness  matrix. 

In  most  cases,  stiffness  matrix  calculations  for  both  BE 
formulations  follow  the  same  main  steps.  However,  the  direct  BEM  is 
better  suited  to  coupled  solution  approaches  because  of  its  simpler 
traction-displacement  relationship.  For  smaller  problems  the 
calculation  of  matrices  comprising  K  require  the  most  numerical  effort. 
For  larger  problems  the  inverse  calculation  of  a  fully  populated  matrix 
dominates  the  numerical  effort.  The  only  known  method  for  reducing  this 
cost  is  to  subdivide  the  homogeneous  BE  region  into  a  number  of  smaller 
regions--trading  the  calculation  of  a  large  inverse  for  the  calculation 
of  several  smaller  ones. 

The  only  matrix  in  the  stiffness  calculation  not  calculated  by 
standard  BE  programs  is  the  shape  function  product  matrix  C.  For 
continuous  BEs  C  is  symmetric.  Additionally  the  locally  based  nature  of 
the  shape  functions  makes  C  very  sparse.  The  calculation  of  C  and  its 
multiplication  are  effectively  performed  at  the  element  level. 


RECOMMENDATIONS 

There  are  several  areas  of  investigation  which  could  improve  our 
fundamental  understanding  of  coupled  solution  techniques  or  improve 
their  numerical  effectiveness.  There  are  also  off-shoot  areas  from  this 
and  previous  investigations  which  are  not  within  the  scope  of  this 
project.  Additionally  there  are  areas  which  must  be  investigated  to 
attain  our  long  term  objective  but  could  require  significant 
implementation  effort.  These  later  areas  are  separately  included  below 
but  would  require  6.2  level  funding.  For  FY87,  the  following  areas  are 
of  interest: 


•  numerical  comparison  of  BEM,  FEM,  st i f fness -BEM,  and  coupled 
so  1 ut ions 

•  numerical  study  of  incompatibility 


between  FEM  and  IBEM  using  the  same  shape  functions 
inherent  to  che  DBEM 

•  theoretical  basis  and  numerical  techniques  for  dynamic  problems 
with  infinite  domains 

•  theoretical  basis  and  numerical  techniques  for  semi-infinite 
domains . 

The  first  two  areas  are  the  recommended  core  for  the  FY87 
investigation.  The  later  two  are  rather  general  and  could  be 
investigated  as  time  permits.  Areas  of  investigation  which  either  are 
6.2  level  efforts  or  require  significant  6.2  support  include: 


•  numerical  study  of  coupled  solution  behavior  for  static  problems 

•  symmetry  considerations  in  formulating  the  BE  equations 

•  further  investigation  of  element  integrations 

special  techniques  for  integrations  over  singularities 
special  techniques  for  the  near  boundary  region 
order  of  integration  calculation 

•  Galerkin  formulation  of  BEMs 

•  development  of  a  DBEM  code  to  support  6.1  efforts 

•  visco-elasticity. 

Two  topics  which  are  "off-shoots"  of  this  and  previous  investigations 


•  adaptive  mesh  refinement  based  on  boundary  condition  violation 

•  BF,  formulated  crack-tip  "elements"  in  FE  hosted  couplings  using 
the  recursive  integration  technique. 

We  appear  to  have  unique  aspects  in  our  approach  to  both  of  these 
prob 1 ems . 
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Appendix  A 


LIST  OF  SYMBOLS 


The  most  general  form  of  a  variable  (i.e.,  with  the  least  amount  of 
supplementary  symbols)  is  defined.  For  example,  u^  represents  a  displace¬ 
ment  vector  but  depending  on  the  notation  it  may  be  continuous  or  discrete 
known  or  unknown.  There  are  some  non-unique  uses  of  Latin  symbols.  For 
these  cases  supplementary  symbols  are  included  as  necessary.  Indicial 
notation  is  only  used  for  vector  and  tensor  quantities.  The  indice  range 
in  this  case  is  two. 

Mathematical  Symbols 

V(x)  Explicitly  defines  the  quantity  as  continuous  as  opposed  to 

discrete. 

V  A  rectangular  or  square  matrix  of  functions  (e.g.,  M(x))  or 
constants  (e.g.,  K) . 

V  A  column  vector. 

V  A  known  value(s);  for  example,  u  is  a  vector  of  known  nodal 
displacements . 

An  unknown  value(s);  for  example,  t  is  a  vector  of  unknown 
nodal  tractions. 

_  T 

Transpose  of  a  matrix  or  vector.  V  is  a  row  vector. 

Inverse  of  a  matrix. 

Latin  Symbols 

A . ^  A  constant  tensor  which  provides  a  zero  displacement 

for  the  fundamental  singular  solution  at  a  reference 

distance  r  .  A.,  =  -6,,  C,[C„  ln(r  )-ll. 

o  lk  ik  1  2  o 

C  Square,  sparse  matrix  which  contains  boundary 

integrations  of  traction  and  displacement  shape 
function  products. 


lO 


Cl,  C2,  C3 ,  C4 


Constants  used  in  defining  the  fundamental  solution. 
Functions  of  the  material  constants. 


e 


k 


<$> 


A  unit  concentrated  force  used  in  the  definition  of 
the  fundamental  solution. 


E 


Young's  modulus. 


E 

E  i  j  k f  -  ’  ^  } 


Square  matrix  relating  nodal  tractions  to  nodal 
displacements,  t  =  E  u 

Kernal  function  in  Somigliana's  2'nd  identity 
relating  boundary  displacements  to  the  stress  field. 


f 


Column  vector  of  generalized  nodal  forces. 


F.  . 

U 


(x,|) 


F 


c..(5,p 


ff 


ijk 


lx, 5) 


Mix) 


Nf  x) 


II 


Fundamental  singular  solution,  Kelvin  solution,  for 
traction . 

A  matrix  of  element  integrations  involving  the 
traction  fundamental  solution.  Both  the  direct  and 
indirect  formulations  consist  of  F  matrices.  The 
two  F  matrices  differ  due  to  the  switch  in  roles  of 
the  field  and  source  point  and  the  indices. 

Fundamental  singular  solution,  Kelvin  solution,  for 
displacement . 

A  matrix  of  element  integrations  involving  the 
displacement  Kelvin  solution.  Due  to  the 
symmetries  in  the  displacement  Kelvin  solution  the 
G  matrices  for  the  direct  and  indirect  formulations 
are  equivalent. 

Kernal  function  in  Somigliana's  2'nd  identity 
relating  boundary  tractions  and  body  forces 
to  the  stress  field. 

Traction  shape  functions  used  with  nodal  traction 
values  to  approximate  a  piecewise  continuous  traction 
d i s  t r ibut ion . 

Displacement  shape  functions  used  with  nodal 
displacement  values  to  approximate  a  continuous 
d i sp 1 acemen t  d is t r i but i on . 

Order  of  the  system  of  equations. 


P .  Artificial  boundary  traction  vector,  analogous  to 

simple-layer  source  or  source  density  of  potential 
proh I ems . 


A -3 


Wx.O 


Greek  Symbols 


Distance  from  the  field  to  the  source  point. 

2 

r  =  y  .  y . . 

1  i 

Actual  boundary  traction  vector. 

Fundamental  singular  solution,  Kelvin  solution,  for 
stress . 

Displacement  vector. 

Position  vector  to  a  point  on  the  boundary  or  within 
the  domain  of  the  problem. 

Vector  from  the  source  to  the  field  point. 

y  .  =  x  .  -  l  .  . 

i  l  l 

Position  vector  used  in  the  integration  of  body 
forces . 


Complete,  finite  boundary  of  the  problem. 


Kronecker  delta  svmbol 


Strain  tensor. 


A  Lame's  constant. 

A  Lame's  constant,  the  shear  modulus. 
Poisson's  ratio. 


Position  vector  to  a  point  on  the  boundary  of 
the  problem. 

Stress  tensor. 

Domain  of  the  problem. 

Vector  of  body  forces. 
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